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Abstract 

We prove approach to thermal equilibrium for the fully Hamiltonian 
dynamics of a dynamical Lorentz gas, by which we mean an ensemble of 
particles moving through a o?-dimensional array of fixed soft scatterers 
that each possess an internal harmonic or anharmonic degree of free- 
dom to which moving particles locally couple. We establish that the 
momentum distribution of the moving particles approaches a Maxwell- 
Boltzmann distribution at a certain temperature T, provided that they 
are initially fast and the scatterers are in a sufficiently energetic but 
otherwise arbitrary stationary state of their free dynamics-they need 
not be in a state of thermal equilibrium. The temperature T to which 
the particles equilibrate obeys a generalized equipartition relation, in 
which the associated thermal energy k^T is equal to an appropriately 
defined average of the scatterers' kinetic energy. In the equilibrated 
state, particle motion is diffusive. 

* Stephan.De-Bievre@math.univ-lillel.fr 
t parris@mst.edu 



1 Introduction 



While it is a well established fact that large systems converge to thermal 
equilibrium starting from a more or less arbitrary initial state, few or no gen- 
eral dynamical mechanisms responsible for this approach to equilibrium have 
been identified. An exception is the case where a small system S with Hamil- 
tonian H$ and only a few degrees of freedom is weakly coupled to a large 
"reservoir" R with Hamiltonian Hr and many degrees of freedom. Under 
these circumstances, the small system will converge from an arbitrary ini- 
tial state po to a state ps ~ exp(—Hs/k^T) of thermal equilibrium provided 
that the reservoir is itself initially in a thermal state pa ~ ex.p{—Hn/k^T). 
This phenomena of return to equilibrium has been proven rigourously in a 
number of (relatively simple) systems {JP981 IBFS00[ IDJ03] . In such sys- 
tems return to equilibrium occurs because, roughly speaking, the smaller 
subsystem acts as a small perturbation of the larger one. Strong stability 
properties of the thermal states p ~ exp(— ////geT) then force the coupled 
system to a joint equilibrium state characterized by the initial temperature 
of the reservoir. Studies of this kind however, do not address the question 
of how the past evolution of the reservoir dynamically led to its degrees of 
freedom being in thermal equilibrium in the first place. 

In this paper we demonstrate approach to equilibrium within a class of 
fully Hamiltonian models that we refer to as dynamical Lorentz gases. In 
these models, each member of an ensemble of particles moves through an 
array of localized independent scatterers, each possessing an internal har- 
monic or anharmonic degree of freedom to which the particle locally couples 
[see (jl.ip ]. The degrees of freedom of the scatterers are initially drawn 
from a stationary but not necessarily thermal state of their uncoupled dy- 
namics. Our main result is that, asymptotically in time, as a result of 
repeated scattering events between a particle and the degrees of freedom of 
the medium, the particle momentum distribution pt(p) will, starting from an 
arbitrary initial distribution po(p) of sufficiently large average mean speed, 
inevitably converge to a Maxwell-Boltzmann distribution characterized by 
a well-defined, non-zero temperature T. This constitutes a true approach 
(rather than return) to equilibrium since this emergence of the Maxwell- 
Boltzmann distribution occurs even when the scattering medium, which 
serves as a reservoir for the particle, is not itself in a state of thermal equi- 
librium to begin with. Our result therefore helps to explain, from a purely 
dynamical point of view, the robustness and ubiquity of the Boltzmann 
factor that characterizes states of thermal equilibrium. 

In our analysis, the effective temperature T that dynamically emerges 
from the interaction between the moving particle and the scatterers leads to 
a generalized equipartition relation [see (|4.1|) ]. in which the thermal energy 
k^T '/2 per degree of freedom of the particle is equal to a suitably-defined 
coupling- weighted average of the kinetic energy of the scatterers. In two 
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particular cases it reduces to the standard (i.e., un- weighted) equiparti- 
tion result: (i) when the coupling of the particle to a scatterer is linear in 
the displacement from mechanical equilibrium of the latter's coordinate (in- 
dependent of the particular stationary distribution of the scatterers), and 
(ii) when the scatterers are themselves initially in thermal equilibrium at 
a well-defined temperature T (independent of the specific coupling between 
the scatterers and the particle). As an example of a situation in which the 
standard equipartition result is not obeyed, we analytically predict, using 
the generalized equipartition relation, and confirm numerically, that when 
the particle is quadratically coupled to "harmonic scatterers" that are each 
initially in their own micro-canonical state of fixed energy, the final thermal 
energy of the particle is exactly one-half the value expected from standard 
equipartition. 

We show finally that in these models the asymptotic mean-squared dis- 
placement of an initially localized ensemble of moving particles grows lin- 
early in time, with a diffusion constant D that depends on the effective 
temperature T through the relation [see (|5.8p ] 

1 5 
D~T U , with - < v < -, 

where the power v depends on the nonlinearity of the coupling and on the 
anharmonicity of the scatterers. A more strongly nonlinear coupling leads to 
a lower value of v, and to slower diffusion. Stronger scatterer anharmonicity, 
on the other hand, leads to a higher value of v, and to faster diffusion. 

Our demonstration of these results depends upon a careful analysis of 
the interaction between the moving particles and the local scatterers, and 
notably of the energy exchange that occurs between them. Based on the 
Hamiltonian nature of the dynamics we show that, for sufficiently fast par- 
ticles, although random fluctuations in the energy change of the particle (per 
scattering event) have a magnitude that is of order the average en- 

ergy change is negative at high energies, and of order ||p|| -2 [see (|3.9p . (|3.7|) 
and (|3.8]) ]. Thus, while random fluctuations tend to induce a diffusive growth 
in the particle's speed to ever higher values, a weaker (but systematic) av- 
erage energy loss acts as a source of dynamical friction |Cha43aj that tends 
to reduce it. The relative strengths of these two processes, as encoded in 
the precise power laws above, are of the exact form required to dynamically 
drive the kinetic energy distribution of the particle asymptotically in time 
to a Boltzmann distribution, with a temperature T determined by the ra- 
tio of the magnitude of the fluctuations to that of the dynamical friction 
[see (|3.24p ]. We stress that, in our analysis, this fluctuation-dissipation-like 
relation emerges naturally from the microscopic Hamiltonian dynamics, and 
acts as a defining property of the temperature of the limiting Boltzmann 
distribution, rather than as an a posteriori property of thermal equilibrium, 
as it does in most treatments (see for example [KTH9lj ). 
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The results presented here constitute a generalization in various ways of 
our previous work [SPB06] in collaboration with A. Silvius. There, diffusive 
behavior was proven for a particle moving through a one-dimensional lattice 
of harmonic scatterers to which the particle was coupled linearly with a 
very particular form factor. Approach to equilibrium, although observed 
numerically, was not addressed in that paper. Our work here benefits from 
the insights gained since in [ABLP10J, where the related but very different 
problem of the motion of fast particles in a random time-dependent potential 
was studied. The potentials considered in that paper, while also generated 
by an array of isolated time-evolving scatterers, are non-reactive (or inert), in 
that they do not respond to the particle's passage. The total energy of such 
a system is not conserved, and the random scattering events experienced 
by the particle cause its kinetic energy to grow slowly but indefinitely, in 
the phenomenon of stochastic acceleration. As discussed above, however, 
when the scatterer degrees of freedom are treated dynamically, as they are 
in the current paper, the Hamiltonian interaction between the scatterers 
and the moving particle provides a way for the particle to dissipate excess 
energy to the medium; this mechanism completely suppresses stochastic 
acceleration and allows the particle to approach thermal equilibrium. In 
fact, in the Hamiltonian models considered here, the frictional component 
of the force cannot be independently made small compared to the strength 
of its fluctuating part. As a result, in such systems stochastic acceleration 
cannot be observed, even on intermediate time-scales, before equilibration 
sets in. 

We now describe more precisely the models that form the focus of our 
analysis. After a dimensional re-scaling of the dynamical variables, masses, 
coupling constants, and time into appropriate dimensionless forms, the particle- 
scatterer system in our model is assumed to obey the following equations of 
motion 

y(t) = -aJ2v(QN(t))Va(y(t) - x N )) 

N \ (1.1) 

MQ N (t) + U'(Q N (t)) = -ari'(Q N (t))a(y(t) - x N ). J 

In these equation^], y(t) S M. d , d > 1 represents the particle position at 
time t, and Qn £ K the displacement of the internal degree of freedom 
associated with the scatterer centered at the fixed point xjv G The 
potential U governing the uncoupled dynamics of the scatterer (i.e., of its 

X A11 variables and constants appearing in this equation of motion can be thought of as 
being dimensionless. They can be obtained from a dimensional model 

= -fil]0(Ojv(T)/^)V<7((y(t)-fi V )/L) 

JV 

&^nr + kU'(Q N (r)/i) = -&r,'(Q N (T)/£)a((y(t) - x N )/L), 
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internal degree of freedom) as well as the coupling function 77 are assumed 
to be smooth functions of their argument. More precisely, we will always 
assume that U and r\ exhibit polynomial growth of the type U(Q) ~ \Q\ r , 
\v(Q)\ ~ \Q\ r > f° r some integers < r' < r, so that in particular U is 
confining. The simplest case, in which rj is a linear function of its argument 
will be referred to as "linear coupling". The locations x^,N G Z d of the 
scattering centers can be chosen either randomly (with uniform density) or 
lying on a regular lattice. The form factor a(-) appearing in the interaction 
terms is assumed to be a rotationally invariant smooth function, bounded 
by one in absolute value and supported in a ball of radius 1/2, so that the 
particle interacts with the scatterer at xjy only when \\y(t) — xn\\ < 1/2. We 
suppose min/v^M \\%N — %m\\ > 1) so that the interaction regions associated 
with different scatterers do not overlap, and further assume that the system 
has a finite horizon, so that the distance over which a particle can freely 
travel without encountering a scatterer, is less than some fixed distance 

> 0, uniformly in time and space and independently of the direction in 
which it moves. Finally, a, a > are coupling constants that we assume 
to be small. A typical numerically computed trajectory y(t) is presented 
in Fig. [TJ for a system of the above type, described more fully in Sec. 01 
in which the scatterers are harmonic oscillators, the coupling is linear, and 
the form factor a is equal to unity inside the circular interaction regions 
indicated, and vanishes everywhere else. 

In the inert models studied in [ABLP10], a = so that the scatterers 
do not respond to the presence of the particle [see 111 this paper we 

are interested in the case where they do, and take a = a. In this case the 
equations (jl.ip are, with p = y and Pjy = MQn, equivalent to the energy 
conserving equations of motion generated by the Hamiltonian 

H = y + y^H scat (Q N ,P N ) + a^r](Q N )a(y - x N ), 

N N 

where 2 

H scat (Q,P) = ^ + U(Q). (1.2) 

Obviously, it is not strictly necessary in our analysis to introduce the sepa- 
rate coupling parameter a. However, as we will see, retention of this extra 

where £ and L are lengths, rh and M masses, k, a, a energies, by introducing 

2 _ k _ ML 2 _ a . _ a 

ml 2 ' ml 2 ' rhui 2 l 2 ' moj 2 l 2 

and performing the change of variables 

t = ojr, y(t) = y(r)/L, Q(t) = Q(r)/£, x N =x N /L, 

which leaves a dimensionless model governed by three independent dimensionless param- 
eters, M, a, and a. 
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Figure 1: The left panel shows an extended part of a typical trajectory, 
starting from the origin, computed for a particle moving through a hexagonal 
lattice of harmonic scatterers, with unit lattice spacing and linear coupling, 
as described more fully in Sec. [H The right panel depicts a close-up view 
of the region indicated by an arrow in the left panel. In the close-up, the 
points where the particle enters and exits each circular interaction region 
are indicated with black dots. 

parameter in the analysis helps to clearly bring out the source of the dy- 
namical friction exerted on the particle by the scatterers, and to aid in its 
identification with the back-reaction of the scatterer. 

It should also be obvious that, in this class of models, evolution in time 
of the internal degree of freedom Q tv (t) of a scatterer does not lead to a 
corresponding change in its location xn, nor in the size or the shape of 
the interaction region experienced by a moving particle that encounters it. 
It does, however, lead to a time-dependent change in the magnitude (and 
possibly the sign) of the interaction energy between the particle and the 
scatterer during the time that the particle traverses the interaction region 
associated with it. When a = and arj = +00, the scatterers form isolated 
and impenetrable hard spherical obstacles; in this limit the model reduces 
to a standard Lorentz gas. In our generalization of such a gas, arj is as- 
sumed to be finite, so the scatterers are soft; a particle with high enough 
energy can then enter the interaction region associated with each one. More- 
over, as noted above, with a = a, the internal degree of freedom of each 
scatterer responds dynamically to the presence of the moving particle in an 
energy-conserving manner, providing a mechanism for the moving particle 
to dissipate excess energy to the scatterer degrees of freedom. It is these 
features, taken together, that motivate our reference to this class of models 
as dynamical Lorentz gases. 
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The rest of the paper is organized as follows. In Sec. [2] we develop a 
coupled random walk description of the evolution of the moving particle's 
momentum and position that forms the foundation upon which the rest of 
our analysis is based. In Sec. [3] we show how this leads to an uncoupled ran- 
dom walk in the particle's energy, which we use to obtain a Fokker-Planck 
equation for the particle's speed distribution function. This allows us to 
prove our main result regarding the approach of the momentum distribu- 
tion to equilibrium. In Sec. 0] we further explore the implications of our 
analytical result for the final particle temperature, and use the result of 
that analysis to make predictions for specific dynamical models. We then 
present numerical results, based on computations of the sort appearing in 
Fig. [H that confirm the basic predictions of our analysis. In Sec. [5] we con- 
sider the particle's motion in position space, show it to be diffusive, and 
calculate the dependence of the diffusion constant on the effective tempera- 
ture T to which the particle equilibrates. We also discuss how the stochastic 
acceleration of the particle, present when a = 0, is suppressed in a reactive 
medium, with a = a. Additional numerical results supporting our calcula- 
tion of the diffusion constant are also presented in that section. Section [B] 
contains a summary, and a discussion of the connection of our results to 
those appearing in early work of Chandresekhar [Cha43a, Cha43b, ICha43c] . 
The Appendix contains details of the perturbative expansions that underly 
our analysis of Sec. [3l 



2 Particle in a field of scatters: a random walk 
description 

We first recall the random walk description developed in [ABLP10] to de- 
scribe a particle moving through an array of random but non-reacting (a = 
0) scatterers, and generalize it to the present situation in which the internal 
degree of freedom of the scatterer responds dynamically to the presence of 
the moving particle (a = a). To that end, we consider a typical trajectory 
generated by the equations of motion (jl.ip for a particle that successively 
encounters scattering centers y n := xn„ at a sequence of instants t n , with in- 
coming momenta p n and impact parameters b n . These quantities are related 
by (see Fig. E} 

1 Pn 

y~ := y (t n ) = x Nn - -e n + b n , e n = --^-r b n ■ e n = 0, ||6 n || < 1/2. 

2 \\Pn\\ 

Note that we view the impact parameter b n as a vector perpendicular to 
the incoming direction. We choose to = and xn =0, with the particle 
starting near a scatterer located at the origin. The initial displacements and 
momenta of the scatterers are assumed to be identically and independently 
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Figure 2: A particle at time t n impinging with momentum p n and impact 
parameter b n on the nth scatterer, centered at the point y n = %N n - 

distributed according to a stationary state 

p(Q, P) = p(H scat (Q, P)) = p(Q, -P) (2.1) 

of the scatterer Hamiltonian ()1.2|) . At each scattering event the momentum 
of the particle undergoes a change Ap n = p n +i — p n . Since no forces act 
on the particle in the region between scattering centers, the momentum 
Pn+i with which it leaves the n-th scatterer is also the one with which it 
impinges on the next. The momentum change Ap n depends on the impact 
parameter b n and on the arrival time t n at the nth scattering center, through 
the displacement Q n := QN„(t n ) and momentum P n := PN n (t n ) of that 
scatterer at the time the particle encounters it. Specifically, we can write 

Ap n = R(p n ,b n ,Q n ,P n ), (2.2) 

where 

r+oo 

R{p, b, Q, P) = -a dr r/(Q(r))V 'a(q(r)) (2.3) 
J o 

in which q(r),Q(r) (i.e., without subscripts) denote the unique solution of 

the single scatterer problem 

q(t)=-ar } (Q(t))V(T(q(t)) i 
MQ(t) = - U'(Q(t)) - aa(q(t)) V '(Q(t)), 
with initial conditions 

g(0) = 6-ip/|H|, P (0)=p, 0(0) = Q and P(0) = P. (2.5) 
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After leaving the support of the nth scatterer, the particle travels a 
distance £ n to the next, which it reaches after a time At n = t n+ i — t n = 
i n /\\Pn+i\\ with impact parameter b n +i, that both depend on the geometry 
of the array of scatterers, and on the dynamics of the nth scattering event 
through the precise point y+ = y n from which the particle leaves the 
n-th scatterer, and through the outgoing direction e n +i. The particle finds 
the next scatterer in the state (Q n+ i, P n +\)- The process then repeats itself. 

Based upon this description of the dynamics, and ignoring the role of 
recollisions, we argue, as in [ABLP10J, that the motion of an ensemble of 
particles through an array of scatterers of this type is well approximated 
by a coupled discrete-time random walk in momentum and position space. 
Each step of the walk corresponds to one scattering event, where the vari- 
ables (Q n ,Pn) that characterize the scatterer, and the variables l n ,b n that 
characterize the approach of the particle onto the scatterer, are drawn from 
the distributions that govern them in the actual system of interest. In this 
manner, starting from a given initial condition (qo,po), the velocity p n +i, 
the location y n +i, and the time t n+ \ of the particle immediately before scat- 
tering event n + 1 are iteratively determined through the relations 



where K n = (b n ,Q n , P n ) and where R(p,n) is defined by (j2.3|) and (|2.4p ~ 
(|2.5p . In this process, the parameters (Q n ,Pn) are assumed to be inde- 
pendently and identically chosen from the assumed stationary distribution 
P (-ffscat (Q, P)), while the b n are independently chosen at each step uni- 
formly from the d — 1 dimensional ball of radius 1/2 perpendicular to p n . 
Without the loss of any essential physics, we have in (|2.6p replaced the ran- 
dom variable t n at each time step with the average distance t* < L* between 
scattering events. As a result ()2.6[) defines a Markovian random walk de- 
termined completely by po and the random process (K n ) n eN- Note that, for 
each n, when p n is defined by (j2.6[) . K n is independent of p n , since the latter 
only depends on values of kj. with k < n. In what follows we write (•) for 
averages over all realizations of the random process K n . In addition, for a 
function / depending on p and k = (6, Q, P), b ■ p = 0, ||6|| < 1/2 we denote 
the average over n by 



p n+ l = p n + R(p n ,K n ) 



tn+1 - *n+ Uphill 
Un+1 = Un+ i*e n+ i 



(2.6) 




(2.7) 



where Cd is the volume of the ball of radius 1/2 in W 



d-i 
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3 Equilibration 



In this section, we show that, for suitable initial conditions described below, 
the distribution of particle speeds for an ensemble of particles evolving ac- 
cording to the random walk (|2.6p asymptotically approaches a Maxwellian 
distribution, and derive expressions for the temperature T to which it equi- 
librates. To this end, we note that the first equation defining the coupled 
random walk in (|2.6p is independent of the remaining two. This allows us 
to temporarily ignore the spatial motion of the particle, and the time, and 
independently focus on the random walk that occurs in the momentum p n 
and/or the kinetic energy E n = ||p n || 2 /2 of the particle as a function of the 
step (i.e., collision) number n. Thus, for example, the energy exchanged 

AE( p ,k) = i (fp + R(p, k)) 2 -P 2 ) 

during a single collision between a particle of fixed incoming momentum p 
and a scatterer, can be written in terms of the solutions to (|2.4p as 

r*+ 

AE = a A(s)a(q(s))ds where A(s) = rj'(Q(s))Q(s), (3.1) 
Jo 

and where i+ is the instant of time at which the particle emerges from the 
interaction region (\\q\\ < 1/2) at the end of the scattering event; note that 
it can be replaced in the integral by +oo. 

In the following we will suppose that the scatterer distribution p has a 
finite mean energy 




#scat(Q, P)p(H SCSLt (Q, P))dQdP 



and that the probability of finding a scatterer with an energy much higher 
than E* is negligibly small. Both the normalized Liouville measure on the 
energy surface H scat (Q, P) = E*, associated with the microcanonical distri- 
bution, and the Boltzmann distribution Z^ 1 exp(— H scat / k^T) , for T large 
enough, satisfy this condition. In addition, we will assume that the parti- 
cles are almost always both energetic and fast. By energetic, we mean that 
the particle has a kinetic energy ||j? n || 2 /2 well above the typical interaction 
potential arj(Q]y n )a(y(t n ) — xjv n ) that it encounters in any scattering event. 
Since rj(Q) ~ Q r , and U(Q) ~ \Q\ r , with r' < r, the typical size of r/(Q) is 
of order El ^ at most, so this first condition can be expressed as 

||p n || 2 » aEl' /r . (3.2) 

In addition, we need the particles to be fast, which means they cross the 
interaction region in a time that is short compared to the typical time 
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over which the scatterer evolves. Since the period te„ of -ff S cat at energy E* 
i _i 

behaves as te„ ~ i£* 2 , this second requirement is equivalent to the relation 

||p n || 2 » El'K (3.3) 

The conditions (|3.2|) and (|3.3|) can of course be imposed on the initial 
distribution of particle momenta, given the distribution p of scatterer en- 
ergies. Provided that, in the distribution to which the particle eventually 
equilibrates, almost all of the particles continue to satisfy these conditions, 
the following analysis will provide an accurate description of the equilibra- 
tion process. Since, as we will show, the particle distribution equilibrates 
to a final temperature T for which ||p|| 2 is of the order k^T ~ P 2 ~ E* 
[see (|3.24p ]. both of these conditions do indeed continue to be satisfied dur- 
ing the particle's approach to equilibrium. 

Since we will focus on fast particles, we expand the function R (p, k) in 
inverse powers of This leads to a corresponding expansion 

for the energy transferred in a single collision, in which the f3^> (k) are ex- 
plicit scalar functions of the incoming collision parameters. We first remark 
that, since t + is of order 1 , it is clear that /?(°) = 0. A straightforward 
computation of fi^ and (3^ 2 \ worked out in the Appendix, leads to the 
following results: 

=ar ] '{Q)^L 



(2) («) ~ 2 [P 2 ri"(Q) - Mrf(Q)U'(Q)] L 



- m (v'(Q)) 2 L 2 (\\b\\), 



(3.5) 



where for k 6 N, and b, e 6 M. d , with e • e = 1, b ■ e = 0, 

L fc (||&||):=^ 1 dAAV(&+(A-~)e). (3.6) 

Note that the right hand side in (|3.6[) is a function of \\b\\ only, as a result 
of the rotational invariance of a. It follows that [see (|2.1|) and (|2.7|) ] 

2 



=0, s i : = = Tj2 (v'P) 2 L h (3-7) 



since p(Q,P) = p(Q,—P). Also, since we assumed that p is stationary for 
the free dynamics of the scatterer generated by -ffscat [see (|2.1|) ]. one readily 
checks that, for any choice of r], 



[P 2 V "(Q) - M7]'(Q)U'(Q)} = 0. 
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Hence 



OLOL 

2M <"' L »» 2 ' 



(3.8) 



Assuming [|p n || is large in the sense of (|3.2p - (|3.3|) . and dropping for the 
moment all higher order terms, we conclude from (|3.4[) that the energy 
change undergone by the particle during the nth scattering event is approx- 
imately given by: 



One recognizes here a dominant fluctuating term (in Hpnll -1 ), independent 
of the value of a, which is of zero average in view of (|3,7p . and a smaller 
subdominant term (in ||pn|| -2 ) which fluctuates about a negative non-zero 
average when a = a. Whereas the random fluctuations of the first term 
have a tendency to increase the particle's speed without bound, the second 
term, while weaker, is systematic and has a tendency to reduce it. It is the 
competition between these two effects when a = a that eventually leads the 
particle to equilibrate, as we show below. If, on the other hand, a = 0, then, 
as shown in |ABLP10] . the particle undergoes a stochastic acceleration, with 
it's speed increasing as \\p n \\ ~ re 1 / 3 . 

Note, therefore, that for this class of systems, we have been able to 
explicitly separate the force acting on the particle as the result of its passage 
through the medium into a frictional and a random part, a notoriously 
difficult problem of statistical mechanics in general (see [KTH91], p. 37). It 
is furthermore clear from the above discussion that the frictional part is due 
to a back-reaction effect: it results from the change induced in the particle's 
motion by the change in the medium's motion, which is itself brought about 
by the passage of the particle. This effect is entirely absent when a = 0. 

We now prove that the two competing effects described above balance out 
so as to drive the particle's momentum distribution precisely to a Maxwell- 
Boltzmann distribution. It is convenient for analyzing the asymptotic be- 
havior to focus on the random walk associated with a new scaled dynamical 
variable 



For that purpose, we first note, using (|3.4p . that 




(3.9) 
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where /% := /3^(re n ). From this last expression we obtain the relations 



lb, 



n+l | 



\\Pn 



3 M 

^Ej^m + o(\\ Pn \n 



\\Pn+l\\ - \\Pn\\ = Yj 

t= 

and consequently 

A|K|| 3 = |K|| 2 A||p ri 



(4) 1 fM) 



{ \\Pnf +1 



\\Pr, 



+ o(ikii- 6 ) 



1 + 



Hp. 



n+l| 



\\P 



n+l | 



1 2 n 



3 R (£) 



i 4 ) + I (f&\ 



3 [Pn' + k /% 



In terms of £ n defined in ()3,10p and the parameters 



+ o(im- 4 ). (3.H) 



p< 



'-n — „ i 
^1 



(3S 4 )-V3^), 



3 I Sf + 2 n i ' 



equation ()3.1ip takes the form 

e„+l = ^ + £n " /in^n V3 + ^nCn 2/3 + 7n£n ' + 0(C 4/3 )- 

From (EH)- (EE) one finds 



(3.12) 



En = 0, 



e 2 n = 1, n n := (j, 



aa 



2(3S 4 )V3m 



(v'Lo) 2 > 0. (3.13) 



(3) 

The rather involved computation of the higher order coefficients j3 n and 
(3^ is performed in the Appendix and leads to the result that [see (|A.10p 
and (IA131) ] 

zaT = 0, 



and that 



_ _ l//?n 4) 

7 " :=7 = 3 Sf 



(3.14) 



where 5f3^ and <5/3^ are defined in (|A.20j) and (|A,14p of the Appendix. 
Using these results (|3.12p can be written 

Cn+l = in + en -7^C 1/3 + 7C 1 + O (C 1/3 ) + 0(C 4/3 ), (3-15) 
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where the notation Oo (||£n|| -1//3 ) means the term is O (||£n|| -1 ^ 3 ) and of zero 
average. Again focusing on particles that over the course of their evolution 
spend an overwhelming amount of time at high speeds, we drop the asymp- 
totically small error terms in this last equation to obtain a one-dimensional 
random walk _ _ 

£n+l=£n + e n TE + 7-- ( 3 - 16 ) 

Comparing this to the random walk for ||p n || 2 m (|3.9|) . one recognizes again 
in the second and third term on the right hand side the effect of the dominant 
fluctuating part of the force and of its subdominant frictional part. The fact 
that, in (|3.16j) . the dominant fluctuating term is independent of the random 
variable £ n itself constitutes an advantage over (|3.9p that will simplify the 
following analysis. It arises from the fact that £ n is proportional to the 
third power of the particle's speed. Again, we stress the all-important strict 
positivity of Jl [see (|3.13p ] that occurs when a = a and that, as we will see, 
is needed to assure that the third term on the right hand side of (|3.16p acts 
as a source of dynamical friction capable of balancing the diffusive growth of 
£ n generated by the independent and identically distributed terms e n of zero 
mean. The crucial role of the much smaller last term in (|3.16|) will become 
evident below. We note that essentially the same equation was obtained 
in [ABLP10] for the non-reactive case for which a = 0, and so Jl = 0. It was 
then proven for that case that £ n ~ y/n, leading to unbounded growth of the 
energy of the particle. This stochastic acceleration is completely suppressed 
when Jl > 0, as we will see. 

The value of 7 in (|3.14j) is obtained as the result of an involved compu- 
tation of the coefficient /?( 4 ) in the Appendix. In particular it is shown there 
that, up to corrections that are small provided a is small and that the bath 
is sufficiently energetic (i.e., that is sufficiently high), 7 depends on the 
dimensionality d of the system through the simple relation 

7 = i(d-2), (3.17) 

which is independent of the model parameters, the interaction rj, and the 
potential U. 

We now study the asymptotic behavior of the random walk f)3. 16j) ex- 
ecuted by the variable £ n for the case in which Jl > 0. For fast particles 
relevant to the present analysis, and small coupling a, the dynamical quan- 
tity £ n oc ||p n || 3 /a 2 will be much greater than unity, and the changes A£ n 
that occur during any given collision will be small compared to the value 
of £ n itself. Thus, over a large number of collisions during which £ n does 
not appreciably change, we may take n to be a quasi-continuous variable, 
in terms of which the random walk (|3.16|) takes the form of a stochastic 
differential equation 

i(n) = e(n) - v' (^n)) , (3.18) 
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where e(n) is the random process satisfying 



(e(n)) = 0, (e(n)e(n')) = S(n 



n 



and where 



'(0 = 2^ ? 2/3 -7lne 



(3.19) 



The corresponding forward Kolmogorov equation for the momentum density 
n) of the particle then reads 



(3.20) 



We point out that (|3.18p can alternatively be interpreted as a Langevin-type 
equation describing the rate of change of the position £ of an overdamped 
particle moving in a confining potential v and subject to a fluctuating ran- 
dom force e(n). In that interpretation the variable n in ()3.18j) plays the 
role of time, and (|3.20p is the Fokker-Planck equation for the probability 
distribution p(£, n) associated with such a Langevin process. It is then clear 
that the solutions of (|3.2U|) converge for large n (i.e., after many collisions) 
to the limiting distribution 



Peq(f) = Peq£o)exp 



2 / d£V(£' 

ho 

AA£ 2 "exp(-37l£ 2/3 ), 



■ Mexp[-2v^)\ 



(3.21) 



where M is a normalization constant. Physically, /5 eq (£)d£ gives, in equi- 
librium, the fraction of collisions for which the final particle speed is 
such that the quantity £ = ^ llPn || 3 /^i has a value lying between £ and 
£ + o?£. Changing variable from £ back to the speed ||p||, one has, in general, 



n) 



! S 1 1 p(^,n), which gives the distribution 



Pcq 



Si 



■Peq(0 



(3.22) 



associated with the fraction of collisions, in equilibrium, in which the final 
particle speed is \\p\\. The average time that a particle emerging from a 
collision with speed \\p\\ remains at that speed after the collision is just 
Thus, the equilibrium probability density p eq ( ||p|| ) , which governs 
the fraction /9 eq (||p[|)dp of time that a particle spends in equilibrium with a 
speed lying between \\p\\ and \\p\\ + d\\p\\, is found from the relation 



Peq 



X 



f/?cq 
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Using (|3.17p . we thus find that the distribution p eq (||p||,t) of particle speeds 
asymptotically approaches a Maxwellian 

Peq(lbll) = M \\ P \\ 1+ ^e W '-M^^ 



2k B T' 

A/q ||pf- 1 exp(-M ; ), (3.23) 



2k B T 

with an effective temperature T given by the expression 



1 r = _(3Si)V3 _ 03d) )2 _ ! r? /2p 2 

2 B 6/Z 



-/3( 2 ) 



2M ^2 



(3.24) 



the meaning of which we shall further analyze in Sec. [H Note the crucial 
role played by the precise value of 7 = \{d — 2) appearing in the last term 
of (|3.16|) ; it produces the correct "density of states" prefactor in (|3,23[) . 
without which the particle speed distribution P eq (||p||) would not strictly be 
Maxwellian. To finally show that the distribution of (vector) momenta p 
approaches the corresponding Maxwellian distribution 

Pcq(p) = 1 —^r exp(-\\p\\ 2 /2k B T), (3.25) 

(2^ B T)2 

it is sufficient to prove that moving particles, as a result of repeated colli- 
sions, asymptotically lose any memory of their initial direction, so that the 
corresponding momentum distribution becomes asymptotically isotropic. 
Such a demonstration is given in Sec. which contains an analysis of the 
rate at which random collision-induced deflections turn the particle, and in 
which that information is then further used to calculate the diffusion con- 
stant that characterizes the growth of the particle's mean-squared displace- 
ment in equilibrium. Before turning to an analysis of the spatial motion, 
however, we consider in Sec. 0] some specific consequences of our analysis of 
the particle's approach to equilibrium, and illustrate some of the predictions 
of our analysis with appropriate numerical calculations. 



4 Generalized equipartition and numerical simu- 
lations 

We begin by noting that the expression for the effective temperature T 
in (|3.24p gives rise to a generalized equipartition relation: for i = 1 . . . d 

1 1 / P 2 \ 1 ri' 2 P 2 

-(ri) T = - kB T = (—^: = — X=-. (4.1) 
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Here (•) T stands for an average with respect to the Maxwell distribution (|3.25|) 
and (•) denotes a weighted average of the kinetic energy of the scatterers 
in the medium, calculated using the weighting function [rj'(Q)] 2 p (Q, P) , 
rather than simply using the stationary scatterer distribution p (Q, P) alone. 
Clearly, in the case in which the coupling is linear, meaning that r/(Q) = Q, 
(|4.ip reduces to the standard equipartition result 



according to which all degrees of freedom of the system have the same mean 
kinetic energy. Note that this result for linear coupling is completely in- 
dependent of the stationary scatterer distribution p : the medium does not 
itself need to be in thermal equilibrium for (j4.2|) to obtain. 

If the coupling is non- linear, so that rf ^ 1, the situation is more com- 
plicated, as we have indicated. Nevertheless, for the special case in which 
the scatterers are themselves in thermal equilibrium at some temperature 
T*, the corresponding canonical distribution Z^. exp(— H Bca t(Q, P)/1cbT^) 
factors into separate distributions for Q and P, and hence, in this case, so 
does the average 



Thus, in this way the dependence on rj disappears, and one again finds 
that the particle momentum distribution converges to a Maxwellian with a 
temperature T = T* given by (|4.2|) . independent of the precise form of the 
coupling. As a result, standard equipartition holds, meaning that asymp- 
totically in time all degrees of freedom of the system have the same average 
kinetic energy, given by ^k^T. 

As an interesting example in which the effective temperature T that is 
reached does not lead to the standard equipartition relation, consider a sit- 
uation in which the scatterers are oscillators with U(Q) = \Q <1 ', each drawn 
from a uniform ensemble with phase space density p(Q, P) ~ 6(H SC a,t(Qi P)~ 
E*) of fixed energy E*, and quadratically coupled to the particle with rj(Q) = 
\Q 2 ■ A simple computation shows that for this case, for alH = 1 . . . d, 



In other words, the mean kinetic energy per degree of freedom of the particle 
equals one half the mean kinetic energy of the scatterer: standard equipar- 
tition does therefore not hold in this case. 

To demonstrate the general features that emerge from the preceding 
analysis, we have performed numerical calculations for particle-scatterer sys- 
tems evolving according to the actual equations of motion (jl.ip , for the two 
cases of linear coupling r](Q) = Q and quadratic coupling r/(Q) = \Q 2 




(4.2) 



if 



,/2p2 — „/2 p2 




(4.3) 
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Figure 3: Approach to equilibrium of the particle speed distribution for a 
particle moving through a lattice of thermally equilibrated harmonic oscil- 
lators with average energy E* = 2, with linear coupling. The solid curve in 
the last panel is a Maxwellian velocity distribution with the predicted final 
thermal energy kT = 2, equal to that of the oscillators. 

discussed above. All numerical calculations were performed with particles 
moving through a two-dimensional hexagonal array of harmonic scatterers 
(U = \Q 2 ), centered at the points xn = N\u + N2V, where u = (a*,0), v = 
(a*/2, \/3a*/2), and N = (N\,N2) G Z 2 . For these calculations the form 
factor was chosen to be a(q) = x(l — 2\\q\\) where x( x ) 1S the usual step 
function, equal to unity for x > and to zero for x < 0, and the lattice 
parameter a* was set equal to 10/9, ensuring a finite horizon. Note that 
with this choice of a form factor, as can be seen in the close-up in Fig. [TJ the 
particle undergoes (easily computed) impulsive changes in its momentum at 
the instants it encounters the edge of an interaction region, but at all other 
times follows straight line trajectories at fixed speed. 

In the numerical results presented in Figs. [3]l6l for a given set of param- 
eters, each particle in an ensemble of 10 4 trajectories was given the same 
initial speed, at a random point on the boundary of the scatterer at the ori- 
gin, with an initial velocity drawn with equal probability from all physically 
possible outward directions. A typical time sequence showing approach to 
equilibrium is presented in Fig. [3] for a system with linear coupling, a = 0.1, 
and oscillators of unit mass (M = 1) initially in thermal equilibrium at a 
temperature such that E* = kT = 2. The histograms in this figure show the 
evolution of an initially sharp distribution of particle speeds (indicated by 
the dashed vertical line) into the predicted 2-d Maxwell speed distribution 
(shown as a smooth solid curve in the last panel of that figure) with a final 
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Figure 4: Approach to equilibrium of the particle speed distribution for a 
particle moving through a lattice of harmonic oscillators each drawn from 
their own microcanonical distribution with average energy E* = 2, and 
linear coupling. The solid curve in the last panel is a Maxwellian distribution 
with the predicted final thermal energy kT = 2. 

effective temperature T* = T = 2, as predicted by the analysis given above. 

A similar sequence is depicted in Fig. 0] for the same system, but with 
the initial state of each oscillator drawn from a uniform distribution on the 
energy surface H scat (Q, P) = E* = 2. The last panel of that figure confirms 
the prediction that for linear coupling the final distribution is the same as 
when the oscillators are thermally distributed. Additional numerical results 
(not shown here) confirm that this limiting distribution is independent of the 
strength of the coupling parameter a as predicted by (|3.24|) . Figures and E] 
show corresponding results for the case of quadratic coupling [rj(Q) = Q 2 /2], 
again with a = 0.1, M = 1, and E* = 2. With quadratic coupling, the 
preceding analysis predicts different final temperatures according to whether 
the oscillators are distributed thermally with average energy E* or uniformly 
on their individual energy surfaces H scat (Q, P) = E* = 2. The numerical 
results confirm this prediction. The limiting temperature of the particles 
that is obtained when the scatterers are uniformly distributed, given in (|4.3|) . 
is one-half the value that emerges when the scatterers are drawn from a 
canonical distribution having the same average kinetic energy. The standard 
equipartition relation does not hold, in this case, since the average kinetic 
energy associated with the particle's degrees of freedom equals one half that 
associated with the oscillator. 
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Figure 5: Approach to equilibrium of the particle speed distribution for a 
particle moving through a lattice of thermally distributed harmonic oscilla- 
tors with average energy E* = 2, with quadratic coupling. The solid curve 
in the last panel is a Maxwellian velocity distribution with a thermal energy 
kT = 2, equal to that of the oscillators. 

5 Diffusion and suppression of stochastic acceler- 
ation 

We now study the motion in position space of an ensemble of particles in 
thermal equilibrium at high temperature, each of which undergoes the ran- 
dom walk (|2.6p under the same hypotheses on the model as in the previous 
section, in particular that a = a. Our analysis closely follows the approach 
of |ABLP10] . where the case a = was treated. At the end of the section we 
compare the behavior of the mean-squared displacement in these two very 
different situations. 

We begin by noting that, since y n +\ = Un + £*e n and e n = an 
understanding of the particle's trajectory requires more than just a knowl- 
edge of how its speed changes in time; one also needs to understand how the 
random changes induced in the particle's momentum by the scatterers cause 
it to turn. To explore this question, we use perturbation theory in (|2.4|) . to 
expand the function 

f^i \\p\r IN 
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Figure 6: Approach to equilibrium of the particle speed distribution for a 
particle moving through a lattice of harmonic oscillators each drawn from 
their own microcanonical distribution with average energy E* = 2, and 
quadratic coupling. The solid curve in the last panel is a Maxwellian velocity 
distribution with a thermal energy kT = 1 predicted by analysis, one-half 
the value found in Figs. [3]|5j 

in inverse powers of In particular, one finds that 

/+oo 
dA Vcj(6 + Ae), so that e • a (1) = 0. (5.2) 
-oo 

Starting from the first equation of (|2.6p and (|5.ip . a simple computation 
[ABLP10] yields the result 

e n +i = e n + S n , (5.3) 

where 

S n = lTli2+0(\\pn\\- 3 ), atf :=a^(e n , Kn ). 
\\Pn\r 

Thus, as suggested in Sec. O an initially fast particle of momentum po under- 
goes small random and independent deflections that cause the unit vectors 
e n to diffuse isotropically over the unit sphere. To determine the number 
of collisions required for the diffusing unit vectors to cover the sphere, we 
compute the value m = m*(j>o) f° r which the correlation function 

m— 1 m— 1 

(||e m - e || 2 ) = ^2 ^2 ($k ■ h') 

k=0 fc'=Q 
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becomes of order 1. Under the assumptions of the random walk (j2.6|) . the 
off-diagonal terms of this last expression clearly have zero average so that 



2\ _ _ H a 



(!) 



2 



eo =m 



\\po\\ 4 



provided that, over n < m collisions, the particle's speed does not apprecia- 
bly change. A more careful computation, carried out in [A BLP10] verifies 
that this is indeed the leading order contribution to the correlation function. 
Hence, we write 

t^cpo) = ML- ^ 

ll a II 

To justify the assumption that the particle's speed does not change apprecia- 
bly during the first m#(po) collisions, recall that, as we have seen, the energy 
distribution of the particle approaches a Boltzmann distribution with a tem- 
perature T given by (j4. If) : since E* is assumed to be large, the same is true 
for ksT and hence the bulk of the particles in the limiting equilibrium dis- 
tribution, which have speeds comparable to the thermal speed vt = ^dksT, 
are both fast and energetic, in the sense of (|3.2p and (j3.3|) . The reasoning 
above then applies if we can argue that during ■= m*(vT) collisions, 
such particles do not appreciably c hange their spee d. To see this, note that 
from (|3.4|) . and from that fact that /3W = = fi( 3 \ we can infer that 



Wp) - ^ + O (W-) - + o ( Wn . (5.5) 

Hence, over a large number n of collisions, an initially fast particle with 
momentum p will, to leading order, decelerate on average at a rate 



dJHI dJHJ gg) 
dn dE W ||p||3- 

On average, therefore, it will take on the order of 

- w - S <5 ' 6) 

collisions for the particle to slow down to speeds comparable to the critical 
speed 

v* = m.ax.{\J aEl ^ r , \l e\ r } 

obtained from (|3.2p - (|3.3p . For particles in equilibrium, which have speeds 
of the order vt, we find, setting tit ■= h*(«t), that 

n T ||« (1) || 2 V 2 n 2 rpl/r 

m T (^(i)) 2 [t]') 2 
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Consequently, at high temperatures the average time it takes a typical parti- 
cle to slow down is much longer than the one it needs to turn, thus justifying 
the assumptions underlying the derivation of (|5.4p . 

We now conclude. In equilibrium, such particles, which will provide 
the dominant contribution to the growth of the mean-squared displacement 
of the ensemble, travel for roughly my collisions before changing direction 
significantly. After this many collisions the direction of motion of such a 
particle will be uncorrelated with its initial direction. The process then 
repeats itself. We therefore expect the large-scale motion of the particle to 
be well approximated by the following random walk: 



In other words, over large enough length scales, the particle motion is es- 
sentially the same as if it traveled on a straight line path for collisions, 
before turning in a random direction and traveling again along a straight 
line over the same distance. Thus, in this last equation, y^ mT denotes the 
position of the particle after a certain number k of these larger excursions. 
With this picture we can then write 



Extrapolating to all values of n, we thus find that (y^) = nm^il, which 
implies (yf) = mxi^vxt, from which we obtain, using (|5.4p . the diffusion 
constant 



The diffusive growth in time of the numerically computed mean-squared 
displacement of an ensemble of particles moving according to the fully 
Hamiltonian dynamics (jl.ip . for the model described in Sec. HI is illustrated 
in Fig. [7] for the case of linear and quadratic coupling. In the numerical 
work presented in that figure, straight lines indicate linear fits to the nu- 
merical data, which are represented by open symbols, and we have again 
taken a = 0.1 and M = 1. The different curves appearing in that figure 
correspond to temperatures such that k^T = 1, 2, 5, 10, 15, 20 and 25. 

We conclude by providing numerical results that show that this oversim- 
plified argument, which neglects the variation in the speed of the particle 
along its trajectory, captures the essential dependence of the diffusion con- 
stant on the model parameters, in particular on the temperature, and on 



the power-laws \rj(Q)\ ~ \Q\ r ' , U(Q) ~ \Q\ r . Noting from that 



U(k+l)m T 



Vkm T + £*rn T e k m T . 





(5.7) 



a (l)||2 _ Q2r' 



,2r'/r 



~ T 



ar'/r 



one finds from (|5.7p that 
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Figure 7: Numerically computed mean squared displacements for an ensem- 
ble of particles moving through a lattice of thermally equilibrated harmonic 
oscillators, for linear and quadratic couplings as indicated, at thermal ener- 
gies k^T = 1, 2, 5, 10, 15, 20, 25. Steeper slopes in each panel correspond to 
higher temperatures. 

1 5 2r' 5 2 , , 
D~T U , -<v = < . (5.8 

2 ~ 2 r - 2 r y ' 

As mentioned in the introduction, this shows that diffusion is enhanced 
by reducing the non-linearity in the coupling (i.e., by lowering r') and by 
more strongly confining the scatterer degree of freedom (increasing r). In 
particular, for quadratic potentials (r = 2), our analysis predicts that D ~ 
T 3 / 2 if the coupling is linear, and that D ~ VT if the coupling is quadratic; 
those predictions are, in fact, as seen in the numerical results presented for 
these two cases in Fig. [HJ which presents the slopes of the linear fits in Fig. [7] 
as a function of the thermal energy at which they were evaluated. 

We conclude this section with a remark on the relation between the in- 
ert Lorentz gas (5 = 0), treated in |ABLP10"] . and the reactive one treated 
here (a = a). Note that n*(p) in (|5.6p behaves as a -1 . This is as expected, 
since it was shown in |ABLP10] that, when a = 0, not only does the parti- 
cle's speed not decrease, it actually increases without bound as \\p n \\ ~ n 1 ^ 6 
as a function of the collision number n. But, as shown in [AB LP10] . for 
initially fast particles of momentum p, this stochastic acceleration begins 
to manifest itself only after n s (p) ~ ||p|| 6 collisions. Thus, the time scale 
n *{p) ~ l|p|| 4 over which dynamical friction manifests itself by slowing down 
fast particles, and which characterizes the equilibration process, is always 
much shorter than is required by the relatively slow process of stochastic 
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Figure 8: Temperature dependence of the diffusion constant D obtained 
from the data in Fig. [7] for the case of linear coupling (r 1 = 1) and quadratic 
coupling (r' = 2). 

acceleration. Thus, in the reactive models considered here, stochastic accel- 
eration is completely suppressed and cannot manifest itself before the system 
equilibrates. This is a result of the fluctuation-dissipation theorem, which 
implies that the frictional component of the force cannot be independently 
made small compared to the strength of its fluctuating part. 

6 Discussion 

In this paper we have identified a dynamical mechanism that induces ap- 
proach to equilibrium in a classical Hamiltonian system. We have demon- 
strated numerically and analytically that as classical particles scatter repeat- 
edly off the internal degrees of freedom of an array of independent oscillators 
with initial conditions drawn from a fixed invariant distribution, their mo- 
mentum distribution approaches a Maxwellian. The effective temperature 
that characterizes this distribution is obtained from a generalized equiparti- 
tion relation that involves details of the nature of the coupling between the 
particle and the scatterers, as well as of the initial invariant distribution of 
the scatterers. Since our analytical results are derived in a high temperature 
and weak coupling regime, they strictly speaking only establish the precise 
behavior of the equilibrium distribution at large values of the momentum. 
The numerical data presented in Figs. OH on the other hand, clearly indi- 
cate that the momentum distribution converges to a Maxwell distribution 
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over the full range of particle speeds, including very low ones. Identifica- 
tion of the dynamical mechanism responsible for this agreement with the 
Maxwell distribution down to the very lowest speeds is an interesting open 
problem. 

We have also, in our analysis, implicitly focused on the behavior of an 
ensemble of particles each of which has an entire medium of scatterers to 
itself, as it were. We have, in particular, not considered effects that might 
arise in a single system possessing a finite density of moving particles that 
would, as they equilibrate by exchanging energy with the scatterers, cause 
the scatterer distribution itself to evolve (presumably) towards a limiting 
distribution. Indeed, one anticipates that the state of an ensemble of scat- 
terers would, upon being subjected to repeated collisions by a sequence of 
independent itinerant particles, similarly evolve towards a thermal distribu- 
tion. The demonstration of such a result would provide an obvious next step 
towards a more complete dynamical demonstration of approach to equilib- 
rium for the dynamical Lorentz gas, treated as a whole. 

The term dynamical friction used in the present work was adopted from 
a series of papers }Cha43aj ICha43bl [Cha43c] by Chandrasekhar, in which 
he studies the escape rate of a tracer star from a galactic cluster, the latter 
being treated as a gas of stars. Chandrasekhar first argues, on general 
grounds, that if there were no dynamical friction, defined as a deceleration 
of the star along its direction of motion, galactic clusters could not exist for 
as long as they do; stars would escape from them too frequently as a result 
of stochastically induced fluctuations in their velocities to values exceeding 
the escape velocity. Assuming the velocity distribution of the stars to be 
Maxwellian, he then shows that a friction term is indeed produced as a 
result of successive two-body scattering events between the tracer star and 
the other stars of the galactic cluster. Finally, he argues that it has a value 
compatible with the hypothesis of a Maxwellian velocity distribution. 

Our analysis and concerns in the present work, although carried out on 
very different models, are clearly very similar in spirit to those of Chan- 
drasekhar's. The main differences between our analysis and that of Chan- 
dresekhar can be identified as follows. Chandrasekhar computes the average 
momentum change parallel to the incoming velocity p of the star, per colli- 
sion, finding Apn ~ — C/||p|| 3 , where C > is a constant that depends on 
the various parameters of the system. The average momentum change along 
the direction of motion is therefore negative and of order ||p|| -3 , and it is 
this mechanism that Chandrasekhar identifies with dynamical friction. 

In the notation of the the current work, this corresponds to a calculation 
of the quantity e-R(p), and indeed for the current models we obtain a 
similar result: e • R(p) = e ■ a( 3 )||p|| -3 + . . . with 

= ~f}V) - I a (i) . a (i) = - — (^)2 If _ laVllVLoP < 0. (6.1) 
2 2 j\ I 2 
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We argue however that strict negativity of e • a( 3 ) does not, by itself, al- 
low one to conclude that the system can maintain, let alone approach, a 
Maxwellian momentum distribution. Indeed, as we saw, for that purpose 
one needs to establish a stronger condition: it is the coefficient (3( 2 \ that in 
our analysis clearly governs the effect of dynamical friction on the particle's 
speed (or energy), which must be strictly negative. 

To emphasize this point we note that, in our models, even when the 
Lorentz gas is inert so that a = 0, the value of e • a( 3 ) is strictly negative, 
as it is in Chandresekhar's analysis. But in the inert Lorentz gas there is no 
dynamical friction, because /3( 2 ) = (see [ABLP10] ); as a result, particles 
in the ensemble do not equilibrate, they undergo an unbounded stochastic 
acceleration of the sort that Chandresekhar sought to avoid. In fact, in our 
models, and probably in Chandresekhar's as well, a negative contribution to 
the average change in the component of momentum along the direction of 
motion arises, at least in part, because particles turn (as we have discussed 
in Sec. [5]), an effect that has nothing to do with friction or the dissipation 
of excess energy. 
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A The high momentum expansion of the energy 
transfer 

In this Appendix we determine the expansion coefficients in (|3.4|) . for 
I = ... 4, as well as their means and variances, which are crucial ingredients 
of the main results derived in the bulk of the paper. We obtain these from 
(|3.ip and a perturbative calculation of the solutions to (|2.4p to sufficiently 
high order. To this end, it is convenient to introduce a notation in which, 
for any k > 

q(s) = q {k) (s) + 0(s k+1 ), 

and similarly for other functions of s. Here, the notation q^)( s ) designates 
any function of s such that q(s) — q(k)( s ) = 0(s fc+1 ); it is, in particular, not 
unique, an observation that will be helpful in the computations below. Also, 
since we always have s < i+, which is of order it follows that any 

term of 0(s k+1 ) is automatically of 0(||p|| fc+1 ). Then, with q = b — |e, we 
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have 

?(3) 00 = q + ps- ar)(Q)VK (s\\p\\, \\b\\)\\p\\- 2 

-a^ n \Q)VK\{» K il,J 1 
where, for b ■ e = 0, we define, [compare f)3.6|> ] . 

L k (jM,\\b\\)= rd\a(b+(\-l)e) dA, 

JO 

K k (fi,\\b\\) = TdAL fc (A, 

JO 



(A.l) 



Hence AE defined in ()3,ip can be written as follows: 

A£ = AE a + AE b + AE c + 0(\\ P \\- 5 ) (A.2) 

where, with A(s) = rf(Q(s))Q(s), 
rt+ 

AE a = a ds A(s)a(q + ps) 
Jo 

AE b = -^SlJ t+ ds A( s )VK ( S \\p\\,\\b\\)-Va(q + ps) 

AEc = " "Miff Io +dS A ( s ) V ^lNMl b ll)- 

We need to expand each of these terms to order ||p|| included in order to 
determine /3^,£ = . . . 4. In fact, as we have noted in Sec. O f3^ = 
because t + is of order Furthermore, since < s < t + and t+ is 

of the order to obtain the contributions j3a , /St , /3c^ to from 

Ai? a , A^fe, and AE C requires that we expand A(s), respectively, to order 3, 
1, and 0. 

We first determine (3 (1 » = (3^,(3^ = (3 { a\ to which AE b and AE C do 
not contribute. To expand AE a , we compute, from the equations of motion 
for Q(t), 

where 

S(s) = -jj j° dr [U\Q m {r)) + a V ' {Q {2) {r))a{q m {r))\ , (A.3) 
and hence 

G<s)00 = Q + ^-if U '(Q)s>-^U»(Q)s* 



a 



;r]'(Q)K (s\ 



M\\p\\ 2 
aP 

- J ^rj"(Q)K 1 (s\\pl\\b\\). (A.4) 
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From this, one finds, with Lq as defined in (|A,1|) . 

Now using the result 

A < \ '(Ci\ P V'^ 11 '^) 

A(i)00 = r?(Q)- S 



+ ^' ,(Q)s -^M- Lo(s||p|l ' l|6|l) ' (A - 5) 



and 



AE a = a f + ds A {1) (s)a(q + ps) + 0(|M|- 3 ), (A.6) 
Jo 



one easily obtains (j3.5|) . 
We now turn to 



0® and ^=0P+ff>+(J<?>, 

the calculation of which is more involved. Fortunately, in order to establish 
(|3.16p only the averages of these quantities, in the sense of (|2.T[) . are needed; 
these are much easier to compute. 

As a first step in this calculation we note that, since to lowest order 
A(0) = r,'(Q)P/M, 



d l P l n'(OY r +oc 

AEc = M*\W I dAV ^( A 'll 6 ll)- V ^ + Ae ) + °(IHI )' 

(A.7) 



Next, using ()A.5p . we easily determine the leading terms of 

AE b = -^P- f t+ dsA {1) (s)VKo(s\\p\\,\\b\\)-Va(q + ps) 



\P\\ Jo 



+0(||p|| 



o(3) o(4) 

h +^+o(bir 5 ) 
\\p\\ 
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where 



p r+oo 

a?ri(Q)r/(Q) — J d\ VK (X, \\b\\) ■ Va(q + Ae) 



tf } = -^ 1 [P 2 V"(Q)-M V '(Q)U'(Q)] 



r+oo 

/ dX XVK {X,\\b\\) -Va(q + Xe) 
Jo 



a 2 a V (QW(Q) 2 
M 

r+oo 

x / dA L (A, \\b\\)VK (X, \\b\\) ■ Va(q + Ae). 
J o 

(A.8) 

Now, to compute (5^ we need to identify the s 2 contribution to A(s). For 
that purpose, we first write (|A.3j) as 

S(s) = S 1 {s) + S 2 (s) + 0(s 5 ), (A.9) 

where 

Si(s) = ~^J S aT[U\Q {2) (T)) + ar,'{Q m (T))a( q + pr)} 



M 
aa 



o 



S2(s) = j^r,(Q) l_ dTr,'(Q {2) (T))VK (\\p\\T,\\b\\)-Va(q + pr) 



aa 



"(W)(Q) dXVK (X, \\b\\) • Va(q + Ae) + 0(s 4 ). 

Jo 



M\\p\\ 3 

Because «Si(s) is of order s and S^s) is of order s 3 , it follows that 
A(a) = ^(Q (2) ( S )) + Sx(8)ri(Q + ^s) + 0(s 3 ). 

Identifying the s 2 term in this expression is now easily done using the rel- 
evant terms of (|A.4p ; one observes that they are linear or cubic in P and 
hence of zero average in any stationary distribution p. In conclusion, since 
/#=0, 



/3(3)=4 3) +^ 3) =0. (A. 10) 

It remains to determine Since we have already computed f3^ and 

[see (|A.8p and (|A.7|) ]. it is sufficient to determine j3^\ To simplify the 
computation, we make the following two observations. First, when a = 
in the equations of motion the particle moves in a time-dependent 

potential arj(Q(t))a(q) where Q(t) solves the one-dimensional equations of 
motion of a particle of mass M moving under the influence of the confining 
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potential U. As a result, Q(t) is periodic in time, and so, therefore, is the 
interaction potential ar](Q(t))a(q). The corresponding results of ABLP10J, 
which were derived under conditions which include this case, can then be 
applied to this situation. In particular, the coefficients /3W(a = 0), for 
i = . . . 4 were computed in [ABLPlO] and found to satisfy the equalitjjf] 

1 



/3(4)(5 = o) = -(d - tyift 1 )) 2 ^ = 0). (A.ll) 



Using ()3.5p this yields 



/3(4)(« = ) = -(d-3)(/3(i))2. (A.12) 



Hence 



/3(4) = /?(4)(a = 0) + <5/3( 4 ) 

= \{d-Z)WW + W ] +W ] (A. 13) 

where we have used the result 5(3^ = which follows from (|A.7p . Moreover, 
5(3^ can be read off from (|A.8|) . i.e., 

^(4) = _ a 2 a V (QW(Q)i J+°° dA ^ m)VKo{ ^ m . Va{q + Ae) . 

(A.14) 

It is therefore sufficient to determine 5f3^\ which turns out to be less diffi- 
cult than to determine fii^ directly. Indeed, and this is our second obser- 
vation, all expansion coefficients that we compute, and in particular /?„ , 
are polynomial in a and a. Thus, to compute 5(3^ it suffices to com- 
pute only those terms that are of non-vanishing order in a. The required 
terms turn out to be of first or second order at the most. Taken together, 
these two observations allow us to compute far fewer terms, and result in 
a considerable simplification that is ultimately brought about by the fact 
that (|A.12|) is neither obvious nor easily established. To prove it by direct 
expansion in the present context would, in fact, necessitate combining terms 
from fi^\a = 0), f3a\ct = 0), and P^\a = 0). A different, more efficient 
approach is used in ABLP10J. 

To compute 5(3^\ we start by rewriting AE a in the form 

AE a = aA(\\p\\,K, a,a)^ + aB(\\p\\, k, a, a) + 0(||p|r 5 ), (A.15) 



2 Here and in what follows, whenever a function / depends on a, a, and possibly on 
other variables such as ||p|| and n, we shall write f(a = 0) for the values of the function 
/ on the hyperplane 3 = 0, and 8f :— f — f(a = 0). 
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where 



A(\\p\\,K,a,a) = jj^J dA Vf (Q(3) °"^ + Ae ) ( A - 16 ) 



and 



B(\\p\\,K,a,a) = -ytj dA ~n' (Q(3) ^ A 



x5 (m) " (<7+Ae) 

«, a, a) + S 2 (||p||, k, a, a) + 0(||p||- 5 ), 

(A.17) 



with 



Si(||p||,/6,a,a) = ^/ dA? ?'( Q (3) (^)) 



x 5 i ( lAi ) °(<I + ^) ( A -!8) 



B 2 (\\p\\,K,a,a) = J ^{ v (r l 'f)(Q) J\\a(q + \e) f\\' 

WK (X',\\b\\)-Va(q + X'e). (A.19) 

Then, using (|A.15p . we find thatH 

K,a,a). (A. 20) 

We first compute 

5A(\\p\\,K,a,a) = A(\\p\\, k, a, a) - A(\\p\\, k, a, 0) 



x«2<3) CX) ^(9 + Ae) + 0(|b||- 5 ). 

Using (|A.4|) this becomes 
<^4((IH|,K,a,d) = 

° r/(Q) [ 1 d\ V "(Q + A.f-)K (X,\\b\\)a( q + Xe) 



M\\pf'^'J Q ' ^ \\p\\M 
a 



4 P V "(Q) 2 1 1 dA JTxCA, [|6[|)(7(g + Ae) + 0(||^|| 5 )- 



M2||p|| 



o 



(A.21) 



3 Here and in what follows, we write for the coefficient of ||p|| k in the expansion 
of a function / in powers of ||p|| 
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One then easily finds that 



6A W = _ « P 
M 2 



V'(QW"(Q) / dA XK (X, \\b\\)a(q + Ae) 



o 



-V (Q) / dXKtiWlbWMq + Xe) 
Jo 



(A.22) 



Next, we compute the contributions from 5B in (|A.15p . The term coming 
from i?2 can be read off immediately from ()A. 19|) : 



(i) 



SB. 



For 5B\, we write 



^(r/(r/') 2 )(Q) f Q dA a(q + Xe)£dX' 

VK (X',\\b\\) ■ Va(q + X'e). (A.23) 



1 f 1 

<5-Bi(l|p||, = ]p[ _/ dA / ( A ) cr 0'? + eA )' 



where 



= ^(Q(3)(A)) >< S ^ + ^(3)(A'« = °)) >< 



J(A) + J(A), 



so that 



<5-E>j 4 ^ (k, a, d) = / 
It follows from (|A.4jl that 
J(A) = 



+ 



<r(g + Ae). (A.24) 



^ (Q(3)( A'« = ° 



x 



P 



M|H|2 
Using the fact that 



V'(Q)K (X, \\b\\) + —.r/'(Q)Kx(X, \\b\\) S, 

+o(\\ P \\- 5 ) 

n v"(QW(Q)K (X, 11611)5! (±) + 0(||p||- 4 ) 



5*i 



1 



A 



WpWJ m\\ p \\ 



[XU'(Q) + a V f (Q)L (X,\\b\\)] +0(||p||- 2 ), 
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this finally yields 

J(3) ( A ) = ^"(QW(Q)Ko(X, \\b\\) [XU'(Q) + ar,'(Q)L (X, \\b\\)] . (A.25) 

To compute J^ 3 \X), we note that, since 5i(A/||p||) is of order we 
have 

i(A) = rf (q (2) (±,a = o)) ^ (A) + o(ibii^). 

Using (|A.4p . one then finds that 

8S 1 (±- a = Q) = ~ [ M dr [{U'(Q {2) (r)) - U'(Q {2) (r,a = 0)) 



a 



'0 



M\\p\\ 



+ar)' '{Q$){T))a{q + prj] 

A 

2? /(Q) /" W dr C/"(Q (2) (r,d = 0))tf (r|b 



|:^ W dr i/(0(a)(r)Mg + pr) + O(||p||- 

( V 'U")(Q) [ X dX' K (X', 
Jo 

V'(Q)L (X, 



M\W 

a 



M\\p\\ 



a P 

t v"(Q)t7 Li{X, 



M\\p\\ 2 1 W7 M~ 
n AU' V ")(Q)L 2 (X, 



2M 2 \\p\\ 3 

AvW')(Q) I dX> K (X',\\b\\)a(q + X'e) 



a 2 



M 2 \\p\\3 

-M^ PV ' (g)L2(A ' l|6|l) + 0(IHr 

^ + # + # + o(| W |*), 



and 



P A 

MM 

2 



1 [Mwm-p^'iQ)]^ 



2M 2 

+0(\\p\\ 
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so that 



J(3)( A ) = rf(Q)6S[ 3) + V "(Q)^X5S[ 2) (X) 



~ [M(U' V ")(Q) ~ P 2 V"'(Q)] A 2 ^ 1} . (A.26) 



One can now compute /3( 4 ) in rtATTSl) by using (|XTi1) . (TOU^ . (1X221) . (lA34j) . 
(|A.23I) . ()A,25j) and (|A.26|) . Obviously, a general condensed expression for 
/?( 4 ) and hence of 7, defined in (|3.14p . is not readily obtainable. However, 
one can easily check that, under the hypotheses of this paper, namely a = a 
small, E* large and \r}{Q)\ ~ \Q\ r ' , U(Q) ~ with < r' < r, 

7 = ^((i-2) + e(a, J Er 1 ). 

To see this, it suffices to examine the various error terms and use the fact 
that [Ql 17 ~ E v J r and P 2 " ~ £*, for all v € N. 

In the case of linear coupling ry(Q) = Q and an even confining potential 
U(—Q) = U(Q), the computation considerably simplifies since one then 

readily sees that 5/3^ = = 5A^ = <5Pi 4) , so that 



0(4) = ±( d _ 3)(/3(i))' + a£B} 4J (A.27) 



Moreover, for this case J(3)(A) = 0, so that 

dA / 

'0 Jo 

One then finds from (|3,14|) that 



5B{ 4) = /'dA J( 3 )(A) = ^C/"(Q) /'dA /V K (X f ,\\b\\). (A.28) 
Jo M Jo Jo 



7 = -(d-2 ) + — V = -W-2) + a — V 
' 6 V ' 3Sf 6 V ; 3£ 2 



f , 2 x 2^ / 1 dA/ A dA^o(AM|b||)c7(g + Ae) 



Since U" ~ |Q| r ~ 2 , it follows that [/" ~ r « P* ~ P 2 so that indeed 
the last term is negligible for large E*. More explicitly still, for the model 
studied numerically in this paper, where d = 2, U(Q) = \Q 2 and a(q) = 
x(l — 2||g , ||), one finds after an explicit computation of the last term above 
that for this particular case 

3-7T 

^ +1 = 1+ 32Sf' (A ' 29) 
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This produces in (|3.23p a density of states factor 



differing by ||p|| 32fc B T from its Maxwell-Boltzmann value. Even at the rel- 
atively small value of the thermal energy (k&T = 2) used in the numerical 
results of Fig. [3l a deviation in the tails of the distribution of this magni- 
tude in the power law multiplying the Boltzmann factor would be difficult 
to observe without considerably better statistics. 
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